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Abstract Kepler’s Equation is solved using an integrative algorithm developed using homotopy 
theory. The solution approach is applicable to both elliptic and hyperbolic forms of Kepler’s 
Equation. The results from the proposed algorithm compare quite favorable with those from 
existing iterative schemes. 


I Introduction 

The problem of relating orbital position and time in the two-body problem requires the solution of Kepler’s 
Equation, i.e., 

M = E - esinE (1) 

where E is the eccentric anomaly, and M is the mean anomaly for the elliptic orbits. The applications of Kepler’s 
Equation can be classified into two categories [1], 

(1) given a probe’s position in the orbit, determine the time when the probe is at that position 
or 

(2) given the time since last periapsis passage, determine a probe’s position in the orbit at that time. 

The first case is rather trivial. By substituting the current eccentric anomaly E into Eq. (1), one can immediately 
determine the mean anomaly A/, and hence t. The second case, however, is more complex than the first one, since 
Kepler’s Equation has no closed-form solution for the eccentric anomaly E. In the literature, most solution approaches 
for E use iterative methods [1-5]. Recently, a non-iterative method was proposed in Ref. [6] 

In this paper, an alternate solution approach based on homotopy theory is presented. In section II, a preliminary 
introduction to homotopy theory is given. The homotopy algorithm for the elliptic form of Kepler’s Equation is 
developed in sections III and the algorithm for the hyperbolic form is discussed in section IV. A summary and 
conclusion is presented in section V. It is not our intent to claim that our approach is superior (or inferior) to the existing 
iterative approaches, but rather to demonstrate that alternative approaches do exist. However, in order to demonstrate 
that our approach is indeed viable, we do make comparisons with some existing iterative methods [1-5]. These 
comparison results are included in sections III and IV. 

II Homotopy Theory 


Mathematical Model 

The homotopy theory introduced in this paper is based primarily on the research work of Zangwill and Garcia 
[7]. Given a system of nonlinear equations to be solved, the homotopy algorithm starts with a simple solution (a 
solution to a similar but easily solved system of linear/non-nonlinear equations) and through a series of integrations, 
reaches the exact solution of the original system of nonlinear equations. 

Let 9k* be an Euclidean n space. Given f(x) : 9k" —►9k"; i.e., x E 9k", we want to find the solution, 
x = [*;, *;] r ,of 
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( 2 ) 


fix) = 0 

The homotopy procedure for solving Eq. (2) can be stated as follows: 

(1 ) Identify an easily solved system of equations, say 

g(x) = 0 (3) 

where g(x)\ 9V — ► Determine the solution jc°of Eq. (3). 

(2) Define a homotopy function //(x,A), such that 

H(x t 0) = g( x) (4) 

H{xA) = /(*) (5) 

(3) Generate a path (from A = 0toA=l) which leads the solution from jc°to x*. 

Homotopy Functions 

Although there exist numerous forms of the homotopy function, the following three are the most commonly used 
forms: 


(1) 

Newton homotopy 

W(JC,0 =/(*)- 0 - A )J{X°) 

(6) 

(2) 

Fix-point homotopy 

H(x,t) = (1 - A)(jc - x°) + tf(x) 

(7) 

(3) 

Linear homotopy 

H(x,t) = A/( x) + (1 - X)g(x) 




= g(x) + X(f{ x) - g{x)) = 0 

(8) 


The homotopy paths start from an arbitrary point, x ° , for either the Newton homotopy function or the fix-point 
homotopy function; hence, they are the most popular forms of the homotopy function currently in use. However, if 
an easily solved system of equations can be identified, and it is very “close” to the original system of equations, the 
linear homotopy function might be a better mathematical model to use. However, the above homotopy forms can be 
converted to each other; for example, /(*) - /t* 0 ) in Eq. (6) and x — x° in Eq. (7) are equivalent to g(x) in Eq. (8). 

Path Following Algorithm and the Homotopy Differential Equation 


The differentiation of any homotopy function with respect to the homotopy variable, A, is 0, since a homotopy 
function is defined to be zero in (x,A); i.e., 

jir •mi -iff 

(9) 


dH _ dH + dHdx = 0 


dX dX dx dX 

Equation (9) is referred to as the Davidenko differential equation in the literature [8j. Manipulation of Eq. (9) yields 


dx _ _ ,dH\- 
dX K dx> 


i m 
dX 


( 10 ) 


which, along with the condition x(0) = jt°, defines an initial value problem for the path from x° to x". Thus, integrating 
Eq. (10) over X from 0 to 1 yields the desired solution of Eq. (2). Let xE.DC. Sfe", and T = { A I 0 :S A ^ 1}. In 

order to guarantee the path existence, the regularity of H or the inversion of the Jacobian matrix, -gj, in Eq. (10) should 
always hold for all (x,A) in H _l , where H _1 is defined as the set of all the solutions of H(x y X)\ i.e., 

H~ x = {(jc,A) e DxA \ H(x y k) = 0} (11) 

Path existence for homotopy methods are guaranteed by Sard’s Theorem for almost all cases. The path existence 
fails only for points which are in a set of Lebesgue measure zero. 

Theorem 1 Sard’s Theorem (Zangwill [7], Ch. 22) 

Let H : C -* % n , where % is the closure of an open set and H be C*. If K £ 1 + max {0 y q - n}, 

then H is regular for almost all £, except for £ on a set of Lebesgue measure zero, where 
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//(•) = e □ 

Corollary 1 Extended Sard’s Theorem (Zangwill [7], Ch. 22) 

Let D C be the closure of an open set, H : D x A -+ % n be G 2 , and F : D -+% n be G 1 , then the following 
three statements are equivalent, 

(1) H will be regular for almost all e. 

(2) H will be regular at X for almost all e. 

(3) F will be regular for almost all e. Q 

For a homotopy method, q = n 4- 1 , hence if k is greater than or equal to 2, then //(; c, t) is regular for almost 
all e. In other words, as long as H(x t t ) C G 2 , Sard’s theorem states that for an arbitrary e, //(*, /) is almost assuredly 
regular. It can be shown that the second derivative of Kepler’s Equation is continuous. This implies that even if the 
integration path bifurcates during the solution of Kepler’s Equation via the homotopy method, a slight perturbation 
of the stating point of the homotopy equation guarantees that a new integration path exists. 

Ill Kepler’s Equation in Elliptic Orbits 

Homotopy Function for Kepler’s Equation: Elliptic Orbits 

Homotopy methods has been used to solve a system of nonlinear equations [7] [8]. To apply the homotopy 
method, first rewrite Kepler’s Equation (Eq. (1)) as 

/(E) ~ E — esinE — M = 0 (12) 

An easily solved equation is chosen as 

g(E) = E - A# 0 = 0 (13) 


M 0 = M + -j . / 5, slnA/ — - 

1 - sin(M 4- e) 4- sin M 

The linear homotopy function is 

H(E,X) = g(E) 4- X(f(E) - g(E)) 

= E - M 0 4* X(M q - M - e sin(E)) = 0 
Differentiating the above equation with respect to the homotopy variable, X yields 

(1 - XecosE )^~ = esinE 
dX 


dE _ esinE 

dX 1 - Xe cos E ^ ^ 

Notice that the denominator in Eq. (16) is nonzero since Aecos E < 1 for all E. Hence the the homotopy integration 
path can never bifurcate. 

An algorithm for the solution of Kepler’s Equation using homotopy methods is as follows: 

Step 1. Identify an appropriate function g(E) (e.g., Eq. (13)) and develop the homotopy function H(E t X) 
( e.g., Eq. (14)). 

Step 2. Form the initial value problem; i.e., differentiate Eq. (14) to obtain Eq. (16) and solve Eq. (13) 
for the initial condition. 

Step 3. Numerically integrate the initial value problem from X = 0 to X — 1. 

Step 4. Stop. 
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Numerical Results 


In order to demonstrate the viability of the homotopy approach, Kepler’s Equation was solved over a range of 
M and e values using the proposed approach and also using the Newton-Raphson method presented in Ref. [1]. The 
mean anomaly, M , was varied from 0 to 2ji in increments of jr/250 and the eccentricity, e, was varied from 0 to 1 
in increments of 0.02. The numerical integration was accomplished using MATLAB’s ode45 function with 
tol — 5 X 10~ 7 , h = hmax = 1 and hmin — 1/20000. The terminal condition for the Newton-Raphson scheme 
was ! / I < 10” 12 , where f(E) = E — e sin is - M. 

The contour plot in Fig. 1 shows the number of integration steps required to solve Kepler’s Equation as a function 
of mean anomaly M and eccentricity e using the homotopy approach. There are two “4+” regions in Fig. 1; the 
maximum number of integration steps for the left “4+” region is 11 and the maximum number of integration steps for 
the right “4+” region is 29. The maximum number of integration steps in the right “4+” occurs when e = 1 ; excluding 
the case of e = 1, the maximum number is 27. The contour plot in Fig 2 shows the number of iterative loops required 
to solve Kepler’s Equation as a function of mean anomaly M and eccentricity e using the Newton-Raphson algorithm. 
In Fig. 2, there also exist two “4+” regions; the maximum number of iterations in the left “4+” region is 10, and the 
maximum number of iterations in the right “4+” region is 48,818. In the right “4+” region, the maximum occurs when 
e = 1 (note: if e = 1 is excluded, then the maximum number of iterations is 957). 

It is important to observe that, for the Newton-Raphson approach, each iteration involves two function 
evaluations, whereas, for the homotopy approach, each integration involves six function evaluations (using ode45.m). 
Improvements in the performance of the homotopy approach can be achieved by an appropriate selection of the 
function g(E). Our preliminary results indicate that the number of integration steps can be reduced significantly (i.e., 
the contour in the Me-plane is predominantly a 1 integration step region). Our experiments here were rather ad hoc 
and we suspect that a more systematic approach may indeed yield an “optimal” choice for g(E). 



Figure 1 Elliptic Case using a Homotopy Algorithm 


310 






Figure 2 Elliptic Case using a Newton-Raphson Algorithm 


IV Kepler’s Equation in Hyperbolic Orbits 
Homotopy Function for Kepler’s Equation in Hyperbolic Orbits 

For hyperbolic orbits, Kepler’s Equation is 

M = esinhF-F (17) 

where Fis the hyperbolic anomaly, and as before, M is the mean anomaly. To apply the homotopy method, Eq. (17) 
is rewritten as 

/(F) = e sinh F - F - M = 0 (18) 

and an easily solved equation is chosen as 

g(F) = esinhF - M = 0 (19) 

Following the algorithm developed in section HI, the linear homotopy function is 

H(F,X) = g(F) + X(f(F) - g(F)) 

= e sinh F-A/ + 2(— F) = 0 (20) 

which yields the following Davidenko’s differential equation. 

J FT* 

( 21 ) 


(ecoshF - k)^r = F 
ak 


Thus, the homotopy differential equation for the hyperbolic case is 

dF _ F 


dX e cosh F - X 

Once again, observe that the denominator of Eq. (22) is nonzero since e cosh F > X for all F. Hence the the homotopy 
integration path can never bifurcate. 
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Numerical Results 


For this case, the homotopy algorithm was evaluated against two Newton-Raphson algorithms; the first 
algorithm was the “conventional” quadratic convergence algorithm and the second was a quartic convergence 
algorithm [3]. The mean anomaly M was varied from 0 to 3 jt in increments of 7r/150and the eccentricity e was varied 
from 1 to 6 in increments of 0.01. The numerical integrations associated with the homotopy algorithm were again 
accomplished using MATLAB’s ode45 function with tol = 5 x 10~ 7 , h = hmax = 1 and hmin = 1/20000 and the 
termination condition for both Newton-Raphson algorithms was I / 1 < 10" 12 , where /(F) = esinhF — F — M . 
The results are shown in Figs. 3-5. 

The contour plot in Fig. 3 shows the number of integration steps required to solve Kepler’s Equation as a function 
of mean anomaly M and eccentricity e using the homotopy algorithm. The maximum number of integration steps 
for the “5+” region is 9. If instead of Eq. (19), the “known-solution” function was chosen as 

g(F) = esinhF - (M + 0.01) = 0, (23) 

the maximum number of integration steps for the “5+” region is reduced to 7. 

The results shown in Figs. 4 and 5 are based on an initial guess determined from 

F„ = ln(^ + 1.8) 

Figure 4 shows the results for the Newton-Raphson algorithm with quadratic convergence. The maximum number 
of iterations in the “6+” region is 21. Figure 5 shows the results for the quartic convergence Newton-Raphson 
algorithm. In Fig. 5, there are two “4+” regions; the maximum number of iterations in the left “4+” region is 11 and 
the maximum iteration number in the right “4+” region is 5. 

Again, it is important to observe that the quartic convergence Newton-Raphson algorithm requires two function 
evaluations for each iteration, whereas the homotopy algorithm requires six function evaluations for each integration. 
Also, the Newton-Raphson algorithm with local quartic convergence [3] requires five function evaluations for each 



Figure 3 Hyperbolic Case using a Homotopy Algorithm 
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Figure 4 Hyperbolic Case using a Newton-Raphson Algorithm 



Figure 5 Hyperbolic Case using a Newton-Raphson Algorithm with Quartic 
Convergence 
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iteration step. Thus, for the Me - plane analyzed in this paper, the homotopy approach appears to have a computational 
advantage over the quartic convergence iterative approach. Again, additional improvements in the performance of 
the homotopy approach can be achieved by proper selection of g(F). 


VI Summary 

In this paper, we present a new algorithm to solve Kepler’s Equation based on homotopy theory. The procedure 
transforms the root finding problem into an initial value problem which can be integrated to yield the desired solution. 
Although no attempts were made to optimize the developed algorithm, it compared quite favorably with some existing 
iterative algorithms. Improvements in the algorithm can be obtained by investigating appropriate “known functions” 
(g(E) or g(F)) and by investigation other numerical integration schemes. 
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